function [err, EPSCresult]=EPSC_Depression()
% Analyze the Depression of EPSPs
% Derived from EPSP_COV
% Expects data to be a series of EPSCs in a train. The intervals and train
% duration are pulled from the stimulus protocol. 
% also note this needs to look at the recovery pulses.....
% Paul B. Manis, Ph.D.
% pmanis@med.unc.edu
% Initial versions: 4/2002, 5/2002.
% 6/2/2002 (for Covariance).
% pulled and modified 9/8/04 for EPSC..
%

%try
sf = getmainselection;
if(sf > 0) 
    pflag = getplotflag;
    QueMessage('EPSC_Depression analysis', 1); % clear the que
    for i = 1:length(sf)
        EPSC_Depression2(sf(i), pflag);
    end;
end;
return;
%catch
watchoff;
QueMessage('Error in EPSC_Depression Analysis routine', 1);
%end;

%-----------------------------------------------------------

function [err, EPSCresult]=EPSC_Depression2(sf, plot_flag, pflag)

global VOLTAGE CURRENT DFILE
global CONTROL
% first generate empty output arrays

err = 0;
EPSCresult=[];

[DFILE, err] = analysis_setup(DFILE, sf); % do default setup.
if(err ~= 0)
    return;
end;

%-----------------------------Begin analysis--------------------------------
QueMessage('EPSC_Depression - Starting');
dat = [];
time = [];

% abstract general information
protocol=deblank(lower(CONTROL(sf).protocol));
rate = DFILE.rate.*DFILE.nr_channel/1000;
[records,pts]=size(CURRENT);

[RL, err] = record_parse(CONTROL(sf).reclist);
frec=min(RL);
lrec=max(RL);

QueMessage('EPSC_Depression - time base');
% compute the time base for plotting (time is [rec, npts] array for EACH record)
time=make_time(DFILE);
tmax=max(max(time));
k = find(DFILE.rate < 1);
DFILE.rate(k) = 50;
RATES = (DFILE.rate .* DFILE.nr_channel)/ 1000; % array of sampling rates, convert to msec

% always compute the time of each stimulus, in seconds from the start of the data
QueMessage('EPSC_Depression - time base 2');
if(DFILE.mode >= 5)
    wz=DFILE.ztime;
    w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
    if(~isempty(w))
        wz(w+1:end)=wz(w+1:end)+24*60*60;
    end;
    zt = (wz-wz(1))/(60);
    cond_baseline = 5; % 5 min baseline
else
    zt = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
    cond_baseline = 3;
end;
TM=zt;

QueMessage('EPSC_Depression - time base 3');
if(DFILE.mode >= 5)
    wz=DFILE.ztime;
    w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
    if(~isempty(w))
        wz(w+1:end)=wz(w+1:end)+24*60*60;
    end;
    ZT = (wz-wz(1))/(60);
    cond_baseline = 5; % 5 min baseline
elseif length(DFILE.ztime) == length(RL)
    ZT = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
    cond_baseline = 3;
else
    ZT = ones(length(RL), 1)*DFILE.cycle;
    cond_baseline = 5;
end;

% access ZT with the record number (1..n) to get the corresponding ztime
% have to get the records first

% Now, get the times when the valves switched (if any...)
% and generate periods with valve 1, 2, 3 or 4.
QueMessage('EPSC_Depression - Valves...');
p=datac('getnote'); % read the current notefile information.
t_sw_valve=[];
n_valve=[];
TL=[];
VL=[];
if(~isempty(p) & length([p.proto]) > 0) % there should be some, but if not, don;'t do much
    % first set of arrays are immediate representations.
    sw_valve=[1 [p(find(diff([p.valve])~=0)+1).frec]];	% valve switch list (records)
    if(length(sw_valve) > 1)
        n_valve=[1 p(find(diff([p.valve])~=0)+1).valve]; % which valve...
        t_sw_valve=ZT(sw_valve); % don't forget offset from start of data..
        % now make long time arrays to match the other arrays.
        for i = frec:lrec
            TL(i)=ZT(i);
            for j=1:length(sw_valve)
                if(i >= sw_valve(j))
                    VL(i)=n_valve(j);
                end;
            end;
        end;
    end;
end;

% Get analysis windows - some automatic computation is done here...
QueMessage('EPSC_Depression - analysis windows');
stim_list=CONTROL(sf).stim_time; % get the array
psp_time=CONTROL(sf).psp_time; % get psp definition array
if(ischar(psp_time))
    psp_time = number_arg(psp_time);
end;
if(ischar(stim_list))
    stim_list= number_arg(stim_list);
end;
dw = CONTROL(sf).deadwin;
dpsp(1) = psp_time(1)-stim_list(1);
dpsp(2) = psp_time(2) - stim_list(1); 
p2=size(psp_time);
npsc = length(stim_list); % get number of pscs to examine

trmp = floor(number_arg(stim_list(1))./RATES - 0.5); % for RMP/ihold determination
psp_time = reshape([stim_list + dpsp(1) + dw, stim_list + dpsp(2)],length(stim_list),2);

for i = 1:npsc % compute time windows for each psc in the train
    t0(i,:) = floor(psp_time(i, 1)./RATES);
    t1(i,:) = floor(psp_time(i, 2)./RATES);
end;

%-----------------------------Begin analysis in earnest --------------------------------
QueMessage('EPSC_Depression - Measuring Ih and rmp');

% measure holding current and "rmp"
% also measure Rin with s1 pulse after ho (ts2)
% and tau by fitting to s1 pulse.
for i = 1:records
    hold_cur(i) = mean(CURRENT(i,1:trmp(i))); % save in array for later usage
    RMP(i) = mean(VOLTAGE(i,1:trmp(i)));
end;
CONTROL(sf).iHold = mean(hold_cur);
CONTROL(sf).Rmp = mean(RMP);
[CURRENT] = FP_artsupp(CURRENT, DFILE, sf);
% Now smooth the Current out a bit
QueMessage('EPSC_Depression - Smoothing');
for i = 1:records
    fsamp = 1000/RATES(i); % get sampling frequency
    fco = 10000;		% cutoff frequency in Hz
    wco = fco/(fsamp/2); % wco of 1 is for half of the sample rate, so set it like this...
    if(wco < 1) % if wco is > 1 then this is not a filter!
        [b, a] = fir_win(8, wco); % fir type filter... seems to work best, with highest order min distortion of dv/dt...
        ismo(i,:) = DigitalFilt(b, a, CURRENT(i,:)')'; % filter all the traces...
    else
        ismo(i,:) = CURRENT(i,:);
    end
end

QueMessage('EPSC_Depression - Finding EPSC peak currents');
for r = 1:records
    ibase(r)=mean(CURRENT(r,1:trmp(r)));
end;
% find EPSC amplitudes for each entry in the train

% for depression series...
% get psc amplitudes for every pulse and record.
for r = 1:records
    ismo(r,:) = ismo(r,:) - ibase(r);
    for i = 1:npsc
        I_psc(i, r) = min(ismo(r,t0(i,r):t1(i,r)))';
    end;
end;

% reduce to a list without failures - must be true for all pulses in consequective records.
r_ok = [];
failthresh = 100;
for r = 1:records-1
    if(all(abs(I_psc(:, r)) > failthresh) & all(abs(I_psc(:, r+1)) > failthresh))
        r_ok = [r_ok r];
    end;
end;



% compute mean current for each stimulus pulse
for i = 1:npsc
    I_m(i) = mean(I_psc(i,r_ok));
    I_var(i) = var(I_psc(i,r_ok));
end
ipscmin = min(min(I_psc));

% and compute the variance and the covariance of the amplitudes for each pair...
vsum = zeros(npsc, 1);
covsum = zeros(npsc, 1);
covar = zeros(npsc, records-1);
var1 = zeros(npsc, 1);
var2 = zeros(npsc, 1);

for i = 1:npsc
    covsum(i) = 0;
    for r = r_ok
        if(i < npsc)
            j = i+1;
        else
            j = i-1;
        end;
        vsum(i) = vsum(i) + 0.5*((I_psc(i,r)-I_psc(i,r+1))^2);
        covar(i,r) = 0.5 * (I_psc(i,r)-I_psc(i,r+1))*(I_psc(j,r)-I_psc(j,r+1));
        covsum(i) = covsum(i) + covar(i,r);
    end;
    var1(i) = mean(covar(i,:).^2)-mean(covar(i,:))^2;
    var2(i) = (3.5/records)*(std(I_psc(i,:))^2)*(std(I_psc(j,:))^2);
end;
rnum = length(r_ok)-1;
vsum = vsum ./ rnum;
covsum = covsum ./ rnum; % mean covariance

Cvq = 0.32; % quantal variability; is CV of mEPSCs - use to correct per Scheuss, eq. 7.
for i = 1:npsc
    switch i
        case 1
            q(i) = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i)/I_m(i+1)));
        case npsc
            q(i) = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i-1)/I_m(i)));
        otherwise
            a = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i)/I_m(i+1)));
            b = (1+Cvq)^2*((vsum(i)/I_m(i)) - (covsum(i-1)/I_m(i)));
            q(i) = (a+b)/2;
    end;
end;
Ncov = - floor((I_m(1)*I_m(2))/covsum(1)); % N is only calculated for the FIRST pair of responses.
Qc = floor(I_m./q); % This is quantal Content for each response

EPSC_Depression.Qc = Qc;
EPSC_Depression.Ncov = Ncov;
EPSC_Depression.Im = I_m;
EPSC_Depression.Ivar = I_var;
EPSC_Depression.npsc = npsc;

CONTROL(sf).EPSC_Dep = EPSC_Depression; % copy the structure over

QueMessage('EPSC_Depression - analysis complete');


%-----------------------------Prepare for plotting--------------------------------
% for plotting, do baselines/stddev.
%----------------------------- plot if figure is set--------------------------------
if(plot_flag >= 0)
    h = findobj('Tag', 'EPSC_Depression'); % check for pre-existing window
    if(isempty(h)) % if none, make one
        h = figure('Tag', 'EPSC_Depression', 'Name', 'EPSP Depression', 'NumberTitle', 'off');
    end
   datac('addwindow', 'EPSC_Depression');
    figure(h); % otherwise, select it
    clf; % always clear the window...
    fsize = 7;
    msize = 3;
    tmax = max(time);    
    
    [m, b, r, p] = linreg(r_ok', -I_psc(1,r_ok)'); % check drift in response.
    yline = m*[1:records]+b;
    drift = 100*((m*records+b)-(m+b))/abs(m+b);
    % plot valve on top 
    % plot Max voltage for EPSP
    subplot('Position',[0.07,0.35,0.45,0.3]);
    sym = {'ks', 'ko', 'k^', 'k+', 'kx', 'rs', 'ro', 'r^', 'r+', 'rx', 'gs', 'go', 'g^', 'g+', 'gx', 'bs', 'bo', 'b^', 'b+', 'rx'};
    for i = 1:npsc
        plot(r_ok, -I_psc(i,r_ok)/1000, sym{mod(i-1,length(sym))+1}, 'Markersize', 2, 'Markerfacecolor', 'k') % data in black
        hold on;
        if(i == 1)
            plot([1:records], yline/1000, 'linestyle', '-', 'color', 'r');
        end;
    end;
    u = get(gca, 'YLim');
    set(gca, 'YLim', [0 u(2)]);
    set(gca, 'FontSize', fsize);
    ylabel('EPSP Amplitude (nA)', 'FontSize', fsize);
    
    set(gca, 'XTickLabelMode', 'Manual');

    % plot the EPSPs...
    subplot('Position', [0.07, 0.7, 0.45, 0.25]);
    plot(time(r_ok,:)', ismo(r_ok,:)'/1000);
    set(gca, 'XLim', [0, max(stim_list)+10]);
    set(gca, 'YLim', [ipscmin*1.1/1000 0.5]);
    xlabel('ms', 'FontSize', fsize);
    ylabel('nA', 'FontSize', fsize);
    set(gca, 'FontSize', fsize);
    
    % plot 1st vs second epsc
    subplot('Position', [0.57, 0.7, 0.18, 0.25]);
    plot(-I_psc(1,r_ok), -I_psc(2,r_ok), 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', 'none');
    xl = get(gca, 'XLim');
    yl = get(gca, 'YLim');
    m1 = mean(-I_psc(1,r_ok));
    m2 = mean(-I_psc(2,r_ok));
    hold on
    plot(xl, [m2, m2], 'linestyle', '--');
    plot([m1, m1], yl, 'linestyle', '--');
    xlabel('EPSC1', 'FontSize', fsize);
    ylabel('EPSC2', 'FontSize', fsize);
    set(gca, 'FontSize', fsize);
    
    % plot mean and std of each epsp in the train
    subplot('Position', [0.8, 0.7, 0.18, 0.25]);
    hp = errorbar([1:npsc], -I_m, sqrt(I_var), 'ko');
    set(hp, 'markersize', 2, 'markerfacecolor', 'k');
    xlabel('EPSC #', 'FontSize', fsize);
    ylabel('EPSC nA', 'FontSize', fsize);
    set(gca, 'FontSize', fsize);
    P = I_m(2)/I_m(1);
    P5 = I_m(5)/I_m(1);
    PLast = I_m(end)/I_m(1);
    NLast = length(I_m);
    
    u = get(gca, 'YLim');
    u(2) = 1000*ceil(u(2)/1000);
    set(gca, 'YLim', [0 u(2)]);
    
    % compute parabolic fit to mean-variance plot as follows:
    % 1. must go through (0,0), therefore,
    %    y = Ax^2+Bx; C = 0.
    % 2. transform to linear:
    %    y/x = Ax + B. (CV vs mean)
    % 3. do linear fit.
    
    [A, B, r, p] =linreg(-I_m/10^3, -I_var./(10^6*I_m/10^3));
    
    % plot of mean vs variance for each EPSP
    subplot('Position', [0.57, 0.4, 0.18, 0.18]);
    plot(-I_m/10^3, I_var/10^6, 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
    xl = get(gca, 'XLim');
    yl = get(gca, 'YLim');
    hold on
    xp = [0:0.05:xl(2)];
    plot(xp, A*xp.^2+B*xp, 'k-');
    
    ylabel('Var (nA^2)', 'FontSize', fsize);
    xlabel('Mean (nA)', 'FontSize', fsize);
    xl = get(gca, 'XLim');
    yl = get(gca, 'YLim');
    set(gca, 'XLim', [0 xl(2)]);
    set(gca, 'YLim', [0 yl(2)]);
    set(gca, 'FontSize', fsize);
    
    % plot of mean vs variance/mean for each EPSP
    subplot('Position', [0.80, 0.4, 0.18, 0.18]);
    plot(-I_m/10^3, -I_var./(10^6*I_m/10^3), 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
    ylabel('Var/Mean (nA)', 'FontSize', fsize);
    xlabel('Mean (nA)', 'FontSize', fsize);
    xl = get(gca, 'XLim');
    yl = get(gca, 'YLim');
    set(gca, 'XLim', [0 xl(2)]);
    set(gca, 'YLim', [0 yl(2)]);
    set(gca, 'FontSize', fsize);
    
    % plot summary of covariances for pairs...
    subplot('Position', [0.07, 0.07, 0.2, 0.2]);
    plot([1:npsc-1], covsum(1:npsc-1)/10^6, 'ks', 'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
    ylabel('Cov (nA^2)', 'FontSize', fsize);
    xlabel('EPSC Pair', 'FontSize', fsize);
    set(gca, 'FontSize', fsize);
    
    subplot('Position', [0.3, 0.07, 0.2, 0.2]);
    plot([1:npsc], -q, 'ks',  'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-');
    yl = get(gca, 'YLim');
    set(gca, 'YLim', [0, yl(2)]);
    ylabel('quantal size', 'FontSize', fsize);
    xlabel('EPSC #', 'FontSize', fsize);
    set(gca, 'FontSize', fsize);
    
    subplot('Position', [0.53, 0.07, 0.2, 0.2]);
    plot([1:npsc], Qc, 'ks',  'Markersize', 2, 'Markerfacecolor', 'k', 'linestyle', '-')
    ylabel('quantal content', 'FontSize', fsize);
    xlabel('EPSC #', 'FontSize', fsize);
    set(gca, 'FontSize', fsize);
    
    % plot textual information abstracted from analysis...
    subplot('Position',[0.74,0.07,0.25,0.25])
    axis([0,1,0,1])
    axis('off')
    text(0,0.9,sprintf('%-12s R[%d:%d] Protocol:%-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Fontsize', 8);
    text(0, 0.8, sprintf('Comment: %s', DFILE.comment), 'Fontsize', 8);
    text(0,0.7,sprintf('Sol:%-12s  Gain:%4.1f LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
    text(0, 0.6, sprintf('Drift: %6.2f%%', drift), 'Fontsize', 8);
    text(0, 0.5, sprintf('quantal size: %6.1f %6.1f %6.1f %6.1f %6.1f', q(1), q(2), q(3), q(4), q(5)), 'Fontsize', 8);
    text(0, 0.4, sprintf('Qc: %d, %d, %d, %d, %d, Sum: %d, NCov: %d', Qc(1), Qc(2), Qc(3), Qc(4), Qc(5), sum(Qc), Ncov), 'Fontsize', 8);
    text(0, 0.3, sprintf('Facil (p2/p1): %6.3f  (p5/p1): %6.3f', P, P5), 'Fontsize', 8);
    text(0, 0.2, sprintf('Depr (plast/p1): %6.3f  Nlast: %d', PLast, NLast), 'Fontsize', 8);
    text(0, 0.1, sprintf('EPSC1 mean: %6.3f  var: %6.4f nA^2', I_m(1)/1000, I_var(1)/10^6), 'Fontsize', 8);
    orient landscape
    drawnow
    
% this should be under control of a plotflag also, but.....
% automatically print this page to the desktop

basefn = DFILE.filename;
rec = DFILE.frec;
pfile =  sprintf('/users/pmanis/desktop/mat_datac-result-figures/%s_%d_epsc_dep.eps', basefn, rec);
print('-depsc2', pfile); % eps, color, level 2 with this filename...

% fig files can be loaded back in and the graphs modified - which is more convenient than
% regenerating the plots.
pfile =  sprintf('/users/pmanis/desktop/mat_datac-result-figures/%s_%d_epsc_dep.fig', basefn, rec);
saveas(gcf, pfile, 'fig'); % as fig for matlab editing later. , color, level 2 with this filename...
fprintf(1, 'Figure saved as: %s\n', pfile);

return;
% control printing and closing of window for automatic runs
    % f = 1 creates plot and leaves it up
    % f = 2 creates plot but closes it when done
    % f = 3 creates plot and prints it and then closes it
    if (plot_flag > 0)
        print -dljet3;
    end
    if plot_flag == 2  
        close
    end
end;

